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Abstract. Hermite methods, as introduced by Goodrich et al. in 1151, combine Hermite 
interpolation and staggered (dual) grids to produce stable high order accurate schemes 
for the solution of hyperbolic PDEs. We introduce three variations of this Hermite 
method which do not involve time evolution on dual grids. Computational evidence 
is presented regarding stability, high order convergence, and dispersion/dissipation 
properties for each new method. Hermite methods may also be coupled to discon¬ 
tinuous Galerkin (DG) methods for additional geometric flexibility Q. An example 
illustrates the simplification of this coupling of this coupling for the Hermite methods. 


1 Introduction 


The computational simulation of wave propagation is central to geophysical applica¬ 
tions, such as seismic imaging and exploration, the modeling of seismic waves induced 
by earthquakes, and problems in structural acoustics. However, the numerical mod¬ 
eling of intermediate frequency waves is known to be challenging for many standard 
low-order methods, requiring a large number of points per wavelength to adequately 
resolve oscillatory behavior. Additionally, the simulation of propagating waves using 
low order methods is typically subject to significant non-physical (numerical) dissipa¬ 
tion and dispersion. High order methods have the advantage of both rapid convergence 
and decreased numerical dissipation |To|[T7| compared to low order methods, especially 
for problems in intermediate frequency wave propagation |14,^. High order methods 
also tend to have a high number of operations per data access, yielding a computational 
structure well-suited to modern computing architectures |19[|21[[^ . 

Hermite methods, as introduced by Goodrich et al. in | |15[ , are high order methods for 
wave propagation which represent the solution using a piecewise polynomial basis by 
collocating the solution and its derivatives on a structured grid. Solution and derivative 
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information at grid nodes is then used to reconstruct and evolve the solutionin time on a 
staggered (dual) grid. Hermite methods are provably stable and high order accurate for 
hyperbolic equations, including problems with varying coefficients. 

Furthermore, though the reconstruction step requires the access of non-local data at 
neighbor nodes, the computation of derivatives then depends only locally on the recon¬ 
structed data at each node. This is advantageous for high order or multi-stage timestep¬ 
ping methods compared to finite difference methods, where neighboring data must be 
accessed each time derivatives are approximated. This structure has also been noted 
to be well-suited for parallel implementations on modern architectures pipi - Hermite 
schemes, which were initially introduced for Cartesian domains, have also been cou¬ 
pled with discontinuous Galerkin (DG) schemes for numerical simulations on complex 
geometries 0. They have also been applied to problems in aeroacoustics 0, electromag¬ 
netics 0, and fluid dynamics 

The Hermite schemes of Goodrich et al. | [T^ are one instance of a broader family of 
methods involving collocation of the solution and its dervatives. Other methods in this 
family include shape-preserving methods and jet schemes which 

use Hermite interpolation in conjunction with semi-Lagrangian techniques to solve ad- 
vective problems. These differ from the Hermite schemes discussed here in terms of the 
characteristic time evolution procedure; however, the analysis and stability of both Her¬ 
mite and jet schemes both rely primarily on properties of Hermite interpolation under 
high order Sobolev seminorms. 

Sections [LT| and present a generalized view of Hermite methods, and motivate new 
one-step Hermite schemes based on variations in the reconstruction procedure. These 
procedures also aim to simplify the implementation and coupling of Hermite and DG 
schemes 0. Section presents numerical experiments which confirm the high order 
convergence and stability of each method for the advection equation in one dimension. 
Section [^extends each method to two space dimensions and includes numerical results 
for the two-dimensional advection and acoustic wave equations. 


1.1 Time evolution 

In this section, we introduce one-dimensional Hermite schemes for the approximation of 
an evolving solution and its derivatives at a collection of points over an interval [a,b) CR. 
Each Hermite scheme presented has a timestep restriction based only on the domain of 
dependence for hyperbolic partial differential equations. For simplicity of presentation, 
we illustrate this using the ID periodic scalar advection equation 

du _ du 
dt ^dx 
u{x,to)=uo{x) 
u{a,t) = u{b,t), 
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where c is a constant advection speed and Uo{x) is a smooth initial condition. Hermite 
methods may be extended in a straightforward manner to non-uniform grids and more 
general systems of equations with variable coefficients (T^ , though these details are omit¬ 
ted for brevity. 

We define first a primary grid O as a collection of K equispaced points 


= Xm = a+mhx, m = 0,...,K—l}, 


where K is the number of grid points on the interval {a,b) and hx = {b — a)/K denotes the 
spacing between the nodes. For periodic domains, we assume that = x^. 

Next, we introduce the interpolation length scale h (distinct from the grid spacing hx)- 
We assume that a smooth function u(x) is well-approximated over some interval {Lm,Rm) 
with size h = Rm — bmhj a degree N expansion around some point Xm- This expansion 
takes the form 

U(X) « Uni{x) = 

;=0 V 



where hx is some spatial length scale, and is typically taken to be some grid spacing in 
practice. The vector u contains Hermite degrees of freedom, which are scaled spatial 
derivatives at Xm 


Uy 


h^x (^) 


j\ dxi 


, y=o,...,N. 


Xm 


For convenience, we express the advection operator applied to Ufn{x) as an expansion 
around x^ 


dUm{x) 
^ dx 


N 


M 



j 


h{ 


;! 


/ dUm{x) 

V dx 



The coefficients Wj are related to Uj through the derivative matrix D 


w = cDu, 



j^i+1 

otherwise. 


0<i,i<N. 


This yields a semi-discrete system for the degrees of freedom Uj 


du(f) 

df 


—cDu(f). 


We approximate u{x,tn+dt) for some timestep df > 0 by solving this semi-discrete 
system. This is achieved in flst l using a temporal Taylor series; assuming an expansion 
centered around time f„. 


N N-j 

Um{x,t) = Y^Y^ijjk 

j=0k=0 


x — x„ 


t-h 

dt 
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and using the Cauchy-Kowalevski relation, the coefficients Vj^ may be shown to satisfy 

.“-*■ 

or more succinctly using matrix-vector notation 

where U(.,fc) refers to the fcth column of the array U. Inserting these coefficients into 
the temporal Taylor expansion and evaluating at time tn+dt yields an update for the 
solution Um{x,tn)- In practice, this is computed using Algorithm]^ ||^ pl09], which may be 
generalized to linear autonomous systems in multiple dimensions. Alternatively, more 
standard time integration techniques may be used to advance the solution forward in 
time. This evolution of the solution from to = tn+dt may be represented by the 


Algorithm 1 Time evolution procedure for ID scalar advection. 
1: procedure Temporal Taylor series evaluation 

2 : W = u” 

3 : for^ = N,N-l,...,0do 

4 : W = W+ cD)u” 

5 : u"+^=W 


application of some update matrix to the degrees of freedom at time f„, such that 




N 




j=0 


X-Xn 


’ V h. 


As noted in | [T^ , for the scalar advection equation, the Taylor expansion in time yields 
an exact evolution of the approximation to uq, so long as the domain of dependence at fo+ 
dt lies within the interval {Lm,Rm). This translates into a degree-independent timestep 
restriction 

h 

cdt /z = min|x^ ^mi¬ 

ff Lfn and Rm are chosen symmetrically around h the interpolation interval reduces to 
{hm/Ryn) — ~h,Xyn +/l). 

The same domain of dependence argument motivates timestep restrictions for "tent¬ 
pitching" space-time finite element methods, which use a coupled discretization in both 
time and space. The stability of the space-time formulation results in a similar causal 
timestep restriction p2| . 
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1.2 Local Hermite interpolation 


The above evolution procedure hinges on a degree N polynomial representation of some 
smooth function u{x) in an interval {Lm,Rm) which accurately approximates the solution 
and its derivatives at the point This is addressed in using degree N Hermite 
interpolation to produce a degree N = 2N+1 reconstruction at specific points Xm- 

The degree N Hermite interpolant of u{x) over (Lm,Rm) (which we refer to as Um{x)) 
is constructed by specifying (N+1) scaled spatial derivatives at the left and right end¬ 
points of the interval (Lm,Rm) 


T hi d)u 


] 


j\ dxi 


«f= 


hi d-' u 

;! dxi 




The resulting expansions around Lm,Rn 


N 




i=0 


x-u 

h. 


N 




;=0 


X-Rn 

hr 


coincide with the value and first N derivatives of u{x) at each point. 

The Hermite interpolant Um{x) is a polynomial of order 2N+1 over which 

is defined by interpolating the (N+1) solution and derivative values at each endpoint. 
Ufn{x) is represented using the following expansion 


^m(^) ■ 


2N+1 

EU; 

j = 0 


X-Xm 

h. 


where uy are values of the solution and 2N+1 scaled derivative values at These 
coefficients may be determined by solving the interpolation problem 

z = 0,...,2N+l. 


d UfYi 

^ 


3 Ujfi 


dx^ 

j dx‘ 

/ 

l-'m 

dx^ 

Rn, 


This results in the system 


■ ■ 


r 1 

. c"" . 

U = 

1 

1 _ 


(1.1) 


where the constraint matrices C^,C^ G ^ ( 2 N+ 2 ) gj^fQj^ce conditions on Um (x) and its 

2N+1 derivatives at the points 


^mn 



i m—1 

;sn(ii-s) 

s=0 


n—m . m—1 

m n (n-s) 

s=0 


n>m 

n<m 

n>m 


n<m. 
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For convenience, we define H 

H 

as the interpolation matrix which maps solution and scaled derivative data at the end¬ 
points of the interval (Lm.Rm) to a degree 2N+1 expansion at the node Variations in 
the construction of H result in different interpolation/reconstruction schemes. 

For the remainder of this work, we refer to the expansion of Um{oc) as a Hermite re¬ 
construction at Xm- This expansion may then be evolved in time a distance of dt using the 
evolution matrix 



2 Hermite methods in one dimension 


Having introduced evolution and interpolation procedures, one-dimensional Hermite 
methods may now be specified. Given some grid O, Hermite approximation spaces are 
associated with each grid. The degree N space K is defined to be 




mr^m+l) 

tn=0 


nc^[a,b]. 


which consists of piecewise polynomials of degree 2N +1 with N globally continuous 
derivatives at each node 6 O. Each u{x) is defined by (N+1) pieces of Hermite 
interpolation data at each node Xm 


TT ^ 

dxi 


; = 0,...,N+1, m = 0,...,K —1. 


For brevity, we refer to the vector of derivatives at the node as U^. 

Finally, in addition to the evolution and interpolation operators and H, respec¬ 
tively, we introduce the restriction operator R G ]r(^+i)x(2^^+i) such that Ryy = Sjj and 
multiplication of a vector by R extracts the first N+1 entries of that vector. 

Hermite methods march forward in time by combining three procedures over one or 
more stages: 

1. Interpolation using points in the primary grid to produce Hermite reconstructions 
of higher degree N centered around points x^. 


2. Evolution of the higher degree Hermite reconstructions at points x^ forward in 
time to t = tn + dt. 


3. Restriction of the solution by truncating the degree of the polynomial expansion at 

Xfn- 
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We emphasize that for each Hermite method, the stable timestep restriction is indepen¬ 
dent of the degree of approximation N. This is due to the fact that, by evolving the 
solution in time using a temporal Taylor series, a single update step of a Hermite method 
may be interpreted as the composition of the exact evolution of piecewise polynomial 
data with a projection in a seminorm which is preserved by the solution In particu¬ 
lar, Hermite interpolation in one space dimension results in the projection onto pi ecew ise 
degree-2N+l polynomial in the seminorm. Thus, as described in Section 1.1 the 


timestep restriction is determined only by the domain of dependence of the equation and 
the interval of the Hermite reconstruction. 

We refer to the collection of of reconstruction points Xm as an auxiliary grid O on 
which the solution is evolved in time. Different schemes use differing combinations of 
interpolation and evolution procedures, which are summarized in TableDifferent Her¬ 
mite schemes may also vary parameters of the interpolation process. For example, the 
Hermite schemes of Goodrich et al. [IS] (referred to henceforth in this paper as Dual 
Hermite schemes) produce a Hermite reconstruction centered between two grid points, 
which are represented over an auxiliary grid consisting of midpoints of the primary grid. 
These reconstructions are then evolved forward in time and truncated. The resulting 
auxiliary grid data may then be used to compute Hermite reconstructions at the orig¬ 
inal primal points, which are then evolved in time and truncated to complete a single 
timestep. 

We introduce here the Virtual Hermite method, which is equivalent to the Dual Her¬ 
mite method with a timestep of size zero on the auxiliary grid. As a result, the inter¬ 
polations to and from the auxiliary grid may be combined, resulting in a step which 
interpolates and reconstructs on the same primary grid. 

The Central and Upwind Hermite methods aim to interpolate and reconstruct on 
the same primary grid through redefinitions of the interpolation operator. The Cen¬ 
tral scheme expands the interval of interpolation, producing a Hermite reconstruction 
at a point x^ information at neighboring nodes while the Upwind Hermite 

scheme uses data from a single neighbor to produce a directional reconstruction. 

While the implementation of the time evolution operator changes from problem to 
problem, the interpolation procedure remains the same irregardless of the equation solved. 
For this reason, we focus on the description of the interpolation procedure for three spe¬ 
cific Hermite methods in the following Sections. Common to each Hermite scheme is a 
timestep restriction which is independent of the polynomial degree N. 

For each one-dimensional Hermite method discussed, the timestep restriction is pre¬ 
sented only for the scalar advection equation. An extension iokxk systems of hyperbolic 
equations 


du _ du 






kxk 


results in timestep restrictions which are identical to the scalar advection equation, except 
that the wavespeed c is replaced by the spectral radius p{A). A similar analysis extends 
these results to variable coefficient problems, and both are described in more detail in 
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Stage 1 

Stage 2 


Interpolate 

Evolve 

Interpolate 

Evolve 

Dual 

Ijn ijn -^fjn 

^m+1 ^ ^m+1/2 


Ijn+l Tjn — 

^m+l/2 ^ '-’m 


Virtual 

jjfi jm V f m 

^m+1 ^ ^m+1/2 


Ijw+1 TJM — 

[ifi/in + dt) 

Central 

jjfi Tin V f m 

[tfi/in + dt) 



Upwind 

Tjn Tjn vT~jn 

[tfi/l^n + dt) 




Table 1: Overview of interpolation and evolution operations over different stages for 
various Hermite method, where refers to the Hermite solution at the nth timestep 
and nzth grid point. Note that the stable time step restriction results in a different dt for 
each method. 



In the following sections, we introduce the Dual, Virtual, Central, and Upwind Her¬ 
mite methods in more detail. 


2.1 Dual Hermite method 


Hermite methods were originally introduced by Goodrich, Hagstrom, and Lorenz in p^ , 
using "primal" and "dual" grids to facilitate Hermite reconstructions at specific points. 
We refer to this specific Hermite method as the Dual Hermite method. The Dual Hermite 
method introduces the auxiliary grid O, which is taken to be a dual or co-volume grid 

0 = = n+(nz +1/2)/Zj, nz = 0,...,K—1} 

such that the nodes of O are staggered a distance oihx/2 between the nodes of O. For 
periodic grids, the dual grid also satisfies x^+i/ 2 +k — ^m+i/ 2 - We associate also an ap¬ 
proximation space to the dual grid 

~ i ® p2]V+l ^Xyn-l/2 ,x^^i/ 2 )\nC \a,h\. 

J 

Each q{x) G is defined by degrees of freedom 


Qm+l/2,/ 


h{ djq 
j\ dxJ 


^m+l/2 


y = 0,...,N+l, 


m = 0,...,X—1. 


In the Dual Hermite method, the interpolation length scale is hx /2, such that the timestep 
restriction on both primary and dual grids is 


cdt<hx/2. 
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The utility of the dual grid comes in considering the Hermite reconstruction at G O, 
which interpolates at the points /2 and Rm — Xm + hx/ 2, such that 

(^Lm/Rfn) (^m 


In other words, nodal data at points on the primary grid is used to produce a Hermite 
reconstruction Um (x) at each point x^ on the dual grid. This also defines the interpolation 
length scale h — hx/'l, implying a time step restriction oicdi<hx/ 2. 

Denoting the vector of nodal data at a point x^±i /2 as Qm±i/ 2 / the interpolation pro¬ 
cedure is 


U^ = H 


Qm-l/2 

Qm+1/2 


where represents a degree 2N+1 expansion at a point on the primary grid. 

Suppose u^{x) G Pq is the solution at timestep n with degrees of freedom The 
Dual Hermite method interpolates the primary grid solution at x^ and x^+i to the dual 
grid and evolves it in time. The degree 2N+1 solution is then truncated, producing 
^ Denoting degrees of freedom for Qm+Y/2' ^ 

Hermite method may be expressed as 


OM+l/2 

^m +\/2 


= RT*H 


U 

U 


n 

m—1 

n ' 

m+1 _ 




In order to update the solution on the primary grid, the process is repeated, except that 
data from the dual grid is transferred to the primary grid before being evolved in time 


u;;+^=RT‘^^H 


Ctn+l/2 

^m-1/2 

0^+l/2 

^m+1/2 


, m = 0,...,X—1. 


The complete process illustrated in Figure 


2.2 Virtual Hermite method 

The Virtual Hermite method is motivated by the fact that two timesteps in the Dual Her¬ 
mite method may be collapsed into a single update step on the primary grid involving 
nodal data at points x^_i,x^,x^+i, bypassing explicit time evolution on the dual grid. 
This requires the formation of the time evolution operator explicitly, which unfortunately 
depends on physical and discretization parameters, and may vary between timesteps for 
nonlinear problems. We propose the Virtual Hermite method to avoid the explicit con¬ 
struction of for the dual grid. 

The Virtual Hermite method is identical to the Dual Hermite method except for time- 
evolution on the dual grid. In the Dual Hermite method, the solution on both the primary 
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O Q) Q) tn O O O 


^m — 1 ^m+1 ^m — 1 ^m+1 

(a) Primary to dual grid (b) Dual to primary grid 

Figure 1: Interpolation procedures to and from primary and dual grids for the Dual Her- 
mite method. Nodes which contribute data to the reconstruction are circled. 


and dual grids is evolved with timestep dt. The Virtual Hermite method skips one evo¬ 
lution step, taking a timestep of dt = 0 on the dual grid. The evolution operator on 
the dual grid then becomes the identity matrix, and the two steps of the Dual Hermite 
method may be collapsed into a single update step 




m—1 
jjn 

Tjn 

ym+lA 


F = RT*H 


RH 

0 

0 

RH 


where and 

The matrix F resembles the co-volume filter analyzed in p3) , which projects the solu- 
tion on a primary grid to and from a staggered dual grid, suppressing spurious gradients 
of the solution on the primal grid. The Virtual Hermite method uses a similar procedure, 
where F maps Hermite data of degree N to Hermite data of degree 2N+1 by transferring 
to and from the dual grid. We expect Virtual Hermite solutions to resemble filtered Dual 
Hermite solutions, which is supported by numerical experiments in Section [3^ 

The Virtual Hermite method obeys the timestep restriction cdt<hx/c. This is the same 
restriction observed for a single step of the Dual Hermite method, since the interpolation 
length scale h = hx/2is the same in both cases. However, by forming F and using a single 
update step on the primary grid, the Virtual Hermite method eliminates the need to ex¬ 
plicitly store and evolve dual grid solutions. We refer to F as an operator with a 3-node 
stencil since degree N data from the three nodes is required to produce 

degree 2N+1 data at as shown in Figure]^ 

We note that it is also possible to produce higher degree reconstructions using the 
same 3-node stencil. One approach is to directly using the degree 2N +1 reconstructions 
at the two dual grid nodes to produce a final reconstruction of degree 4N+3 at instead 
of truncating the dual grid reconstructions. Another option is to directly interpolating the 
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o 


^m — 1 



O ^n+1 


^m+1 


tr 


^m —1/2 ^m+1/2 


(a) Virtual dual nodes 


O O O ^n+1 



© 

© 

© tn 

^m — 1 


^m+1 


(b) Collapsed dual nodes 



Figure 2: Virtual Hermite interpolation procedure, which transfers to and from an aux¬ 
iliary grid (left) or as a single reconstruction step involving a three-node stencil (right). 
Nodes which contribute data to the reconstruction are circled. 


3(N+1) solution and derivative values at each node to produce a degree 3N+2 recon¬ 
struction at Xm- In both cases, the higher degree reconstruction may then be evolved in 
time using a higher order scheme; however, we do not observe significantly improved 
convergence rates under such a procedure, and for some values of N, higher degree re¬ 
constructions results in an unstable scheme. Section |3^ describes an alternative way to 
determine a higher degree reconstruction based on optimization of discrete dispersion 
and dissipation relations p^ . 


2.3 Central Hermite method 

While the Virtual Hermite method removes the need to update the solution in time on 
the dual grid, the Central Hermite method sidesteps the use of an dual grid altogether by 
defining the Hermite reconstruction at through interpolation at neighboring points 

{Lm-f^m) — {^m ~hx/Xfn -\-hx) — (x^_i,X^+i). 

As a result, the interpolation length scale ish = hx and the timestep restriction is 

cdt<hx. 


Then, nodal data at Xm is constructed using data from nodes at In other 

words, Hermite interpolation at a node on the primary grid is performed using its two 
neighbors and then evolved in time, resulting in an update step 


U^ = RT^^H 


Urn—1 
Um+1 
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O ® O ^n+1 



H 


© 

• 

© tn 

^m — 1 


^m+1 



Figure 3: Central Hermite interpolation procedure, where a Hermite polynomial is con¬ 
structed centered at Xj using nodal information from and X/_i. Nodes which con¬ 
tribute data to the reconstruction are circled. 


The Central Hermite method may thus also be interpreted as two decoupled Dual Her¬ 
mite methods on grids of size 2hx. The timestep restriction and numerical results in 
Section]^ also confirm this interpretation. 


We note that by increasing the size of interpolation interval h, the resulting timestep 
restriction increases independently of the grid spacing /Zj. However, doing so also de¬ 
creases the quality of the interpolation procedure, and Section 3.2 describes deleterious 
effects on the error and spectra of the resulting method. 

The Central Hermite interpolation operator results in the 2-node stencil of Figure 
since nodal information from x^-i and x^+i is required to construct information at a 
Xm. Computationally, a smaller stencil results in fewer memory accesses for the recon¬ 
struction. While the difference between the Central Hermite and Virtual Hermite stencils 
in one space dimension is small, the difference becomes more pronounced in multiple 
dimensions. For a degree N Hermite method in d dimensions, each node in the stencil 
requires (N+1)^ accesses, and a Central Hermite stencil contains 4 nodes in 2D, and 8 
nodes in 3D, while the Virtual Hermite stencil contains 9 nodes in 2D and 27 nodes in 3D. 


2.4 Upwind Hermite methods 

Each Hermite method presented has utilized a centered stencil, where solution values 
and derivatives are interpolated in a symmetric fashion around the reconstruction point. 
The Upwind Hermite method constructs instead a directional or one-sided Hermite re¬ 
construction. This concept was used in to enforce boundary conditions, though the 
use of such reconstructions may also take advantage of the directional nature of hyper¬ 
bolic equations ||^[^. 

As shown in Figure]^ an Upwind Hermite reconstruction Um (x) may be defined at 
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by interpolating solution and derivative values at the endpoints of the interval 


3 Ujfi 



d^Ufn 

^ Ujfi 

dx’ 

X 1 

A.m—1 

/ 

Xfn—1 

dx‘ 

X 


i = 0,...,2N+l. 


The solution to this problem results in a degree 2N+1 reconstruction at Xm- However, 
since the solution and derivative values at Xm are interpolated, the first (N+1) coeffi¬ 
cients at Xm remain unchanged, and only the remaining (N+1) coefficients need to be 
computed. This may be done by multiplying the derivative values at Xm-i,Xm by the last 
(N+l) rows of the interpolation matrix H. A downwind interpolation operator may be 
defined in a similar manner using information at (xm,Xm+i) to produce a reconstruction 

3.t Xjfi. 



(a) One-dimensional reconstruction 


Figure 4; Interpolation procedure in one space dimension, where a Hermite polynomial 
is constructed centered at Xm using additional nodal information from and Xm-i- Nodes 
which contribute data to the reconstruction are circled. 


For the advection equation specifically, the domain of influence is biased, such that the 
solution at time tn+dt at the point Xm depends only on Xm—cdt at time Redefining the 
width of the interpolation interval for the Upwind Hermite method ash = Xm — Xm-i=hx, 
the same domain of dependence arguments used previously imply that the method is 
stable if cdt <hx. The timestep restriction of the Upwind Hermite method then matches 
that of the Central Hermite method, with the caveat that this result is specific to scalar 
advection equations and the sign of c. For example, for c>0, a downwind reconstruction 
would be unstable due to the fact that the interval does not contain the domain 

of dependence for Xm at any time greater than We also note that, by similar arguments 
made in the time evolution of the Upwind Hermite method by temporal Taylor 
series is also exact. 

The Upwind Hermite method also requires special treatment when directionality is 
not readily apparent, such as for systems of hyperbolic equations. In one space dimen¬ 
sion, the procedure may be adapted to reconstruct upwind and downwind characteristic 
variables, similar to the approach used in WENO reconstructions |25 However, the 
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effectiveness of the characteristic approach does not appear to extend to all systems of 
equations in higher dimensions, as discussed in Section]^ 


3 Numerical experiments in ID 

To compare the performance of the new Hermite methods, we examine convergence rates 
and qualitative behavior for the Virtual. Central, and Upwind Hermite methods. Numer¬ 
ical results are shown for the periodic constant-coefficient scalar advection equation on 
the interval [—1,1], using the Taylor expansion discussed in Section p~d] to evolve in time. 

We introduce also a CFL constant C > 0 such that dt = Ch/c, where h is the size of 
the interpolation interval (for Virtual Hermite methods, h = hx/2, while for Central and 
Upwind Hermite methods, h = hx). C 1 sets the timestep as large as possible based 
on the timestep restriction for each method, while C <C 1 results in more timesteps than 
necessary as implied by stability. Since a filter-like step (the Hermite reconstruction) is 
applied at each timestep, small values of the CFL constant C (i.e. smaller timesteps than 
strictly necessary) may result in stronger filtering than necessary and larger errors. 


3.1 Convergence rates 

We report convergence rates for the one-dimensional scalar advection equation with 
speed a = l and solution 

u{x,t)=sm{n{x — t)). 

We vary the CFL constant between C = .1, C = .5, and C —.9 and calculate errors at time 
T = 10. For both methods, the error is smaller the closer C is to 1 as shown in Figure]^ 
Convergence rates are reported in Table and except for the lowest order case N — 1 
and smallest CFL constant C = .1, optimal rates of convergence were observed 

for all methods. Additionally, at higher orders of approximation, errors for the Virtual 
and Upwind Hermite method are very similar to those of the Dual Hermite method. For 
the Central Hermite method, the error is roughly a factor of 2^ greater than that of the 
other methods. 



C = .l 

C = .5 

n 

II 

N 

1 

2 

3 

1 

2 

3 

1 

2 

3 

Dual 

2.72 

4.99 

7.02 

2.93 

5.0 

6.98 

3.02 

5.02 

7.01 

Virtual 

2.67 

5.0 

7.00 

2.96 

5.0 

7.02 

2.99 

5.01 

7.07 

Central 

1.71 

4.94 

7.06 

2.62 

4.98 

6.92 

2.92 

4.98 

6.99 

Upwind 

2.94 

4.99 

6.98 

2.96 

5.0 

7.02 

3.02 

5.03 

7.03 


Table 2: rates of convergence for Hermite methods for the one-dimensional advection 
equation with a smooth sinusoidal solution. 
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(c)C = .9 


Figure 5: Convergence of errors for Hermite methods for the one-dimensional advec- 
tion equation with a smooth sine solution. 


The growth of error in time is also examined for each Hermite method. Error esti¬ 
mates for discretizations of transient hyperbolic problems typically contain two terms 
which characterize spatially-dependent and time-dependent errors, respectively. Stan¬ 
dard bounds are of the form 

e(T)<(Ci + C2T)/z''(^) 

where e(T) is some measure of error at time T, and r(N) is some rate of convergence 
depending on the degree of approximation |8}[T7[. We confirm the linear growth of error 


in time for all Hermite methods in Figure p 
which is discussed in Section [3^ The growt 


with the exception of the case when C = 1, 
1 of error for the Upwind and Dual Hermite 


methods is very similar. While the Central Hermite method develops larger errors than 
the Virtual Hermite method, the long-time rate of growth is identical for each value of 
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C for both methods. Moreover, the Central Hermite method with K = 32 results in time- 
dependent errors very similar to the Virtual Hermite method for K = 16, indicating a 
strong dependence of the error on the size of the interpolation interval. 




(a) Dual Hermite (b) Virtual Hermite 




(c) Central Hermite (d) Upwind Hermite 

Figure 6 : Growth of error in time for various Hermite schemes with N = 3, K = 16 in 
one dimension. The advection equation is solved up to time T — 10 with a smooth sine 
solution, and the error is computed at each timestep. 

Growth of the error in time is often described in terms of dispersive and/or dissipa¬ 
tion mechanisms intrinsic to numerical methods | [^ . This may be illustrated by advect- 
ing an under-resolved function over several periods; the effect of numerical diffusion will 
be to smooth the profile out as time increases. Figures shows advection of the periodic 
Gaussian pulse initial condition ^- 4 sm( 7 rx )2 5 periods for orders of approximation 

N = 1,2,3 and a grid of 8 nodes. 

As expected, all methods display diffusive behavior at low orders of approximation 
and low values of C, which is improved as N and C increase. The Dual and Upwind 
Hermite methods appear to be the least diffusive, though the difference between each 
method is small at higher orders of approximation. On a coarse K = 8 mesh, the Central 



































17 



(a) Dual,C = .l 



(b) Dual,C = .5 



(c) Dual,C = .9 





(j) Upwind, C = .l (k) Upwind, C = .5 (1) Upwind, C = .9 


Figure 7: Advection of a periodic Gaussian pulse by Hermite schemes with various CFL 
constants C on a grid of K = 8 cells. 
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Hermite scheme behaves particularly poorly, displaying both spurious oscillations and 
significant numerical diffusion for all C. 

Increasing to a finer mesh 16, the Central Hermite scheme behaves comparably to 
the Virtual Hermite scheme. Qualitatively, the behavior of the Central Hermite scheme 
for K=16 resembles that of the Dual Hermite scheme for K=8. The errors for advection 
of a Gaussian, while not identical, are very close — for N=2, the Central Hermite scheme 
with K = 16 results in an error of .0841903, while the Dual Hermite scheme with K = 8 
results in an error of .0839913. This is expected since, for Central with K = 16 and Dual 
with K = 8, the timestep restrictions and interpolation intervals are identical. 

3.2 Spectra and dispersion/dissipation relations 

Numerical experiments confirm the high order convergence of each method; however, 
the qualitative behavior of each method in convecting an under-resolved solution varies 
significantly. We seek to further analyze this behavior by computing the spectra of the 
update matrix and dispersion/dissipation relations for each Hermite method. 

We define the update operator S such that, for solution degrees of freedom at time 
G, the application of S evolves the solution at time tn+dt 

Both the Dual and Central Hermite method march forward by dt = Chx/c over a single 
timestep (the Dual Hermite method defines dt = Chx/ (2c), but takes timesteps on both 
primary and dual grids). Since the Virtual Hermite method takes the timestep to be 
dt = Chx/ (2c), we analyze for the Virtual grid (corresponding to taking two timesteps 
instead of one) in order to normalize how far in time the update operator evolves the 
solution. 

For the Dual, Central, and Upwind Hermite methods, setting C = 1 exactly (with 
respect to machine precision) results in exact evolution of the solution, though this is 
unique to constant coefficient equations. As a consequence, the update operator A be¬ 
comes exactly equal to a circulant shift matrix. As a result, the order of convergence of 
each of these methods with C = 1 increases to 0(/7^^+^), and coincides with Hermite in¬ 
terpolation estimates given in Lemma 3.1]. This same exact evolution property does 
not hold for the Virtual Hermite method. 

The spectra in Figure [^suggest that the Virtual and Dual Hermite methods should be¬ 
have similarly, as the eigenvalues of A are distributed similarly for both methods. Eigen¬ 
values which lie on the unit circle are typically of the form and are related to the 
non-dissipative propagation of modes of the form For example, for C = .1, the 

eigenvalues fall closest to the unit circle around the point (1,0), corresponding to the 
non-dissipative propagation of modes with small co (low frequency modes). The remain¬ 
ing spectra lie within the unit circle, indicating dissipation of under-resolved modes. In 
constrast, the spectra for the Central Hermite method clusters not only around (1,0) but 
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(j) Upwind, C = .l (k) Upwind, C = .5 (1) Upwind, C = .9 


Figure 8: Spectra of the update matrix for each Hermite scheme at various CFL constants 
C. The order of approximation and grid size are fixed to be N=3 and K=16, respectively. 
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also around ( — 1,0), suggesting that under-resolved high frequency modes may be prop¬ 
agated without dissipation. These spurious modes may explain the behavior of the Cen¬ 
tral Hermite method observed in Figure where propagation of a Gaussian on a coarse 
grid resulted in "spurious" oscillatory behavior which remained over several periods of 
advection. 

A study of the dispersion and dissipation error for the Dual Hermite method was 
reported in using a modified equation and Bloch wave analysis in one dimension, 
which we adapt and apply to the Hermite methods introduced in this work. Dispersion 
and dissipation properties of Hermite methods depend mainly on the properties of one¬ 
dimensional Hermite interpolation matrix H 




where act on Hermite data associated with a given node and it's left/right 

neighbors to produce a reconstruction. 

For a periodic grid, S G x^(^+i) is a block tridiagonal matrix 


Sl Sr ... Sl 
Sl Sc Sr 


S = 


Sl Sc 


Sr 


Sl Sc 


where Sl = THl, and similarly for Sq^Sr. For the Central and Upwind Hermite methods. 
Sc and Sr are zero, respectively. 

We perform a fully discrete Bloch analysis to examine dispersive and dissipative 
properties of each Hermite method. This is done by representing the wave solution 
e^k{x-ct) Hermite basis, and noting that the solution is shifted in both space and 

time by scaling with a complex exponential 

u{x,t+dt) — u{x,t)e~^^^^\ u{x+h,t) — u{x,t)e^^^. 


Assuming a uniform grid spacing h, the discrete evolution of the interpolated exact solu¬ 
tion at a node x^ from time tn to tn + dt is then given by 

Since the timestep restrictions for the Dual and Virtual Hermite methods are di — Ch/ (2c) 
as opposed to di — Ch/ c, the dispersion relations are measured over two timesteps. For 
the Dual Hermite method, this implies that S captures the evolution of the solution from 
the primal to dual grid, then back to the primal grid. For the Virtual Hermite method, 
this requires taking two timesteps and substituting for S the matrix S^, which is block 
pentadiagonal. 
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Figure 9: Dispersive and dissipative error 
scheme, with N = 1,2,3. The computed errors 


10 ° 



kh 


(b)C = .9 

I Ah — I /1 I £qj. Hermite 

e observed to behave as O {kh/ 


An eigenvalue problem may be solved for the discrete Floquet multiplier A^, whose 
real and imaginary parts correspond to numerical dispersion and dissipation, respec¬ 
tively. For each Hermite method, we measure the relative error between the discrete and 
the true Floquet multiplier 

over a single timestep as a function of kh, the order of approximation N and the CFL 
constant C = .1,.5,.9. 

Figure]^ shows the error E^h over a range of kh. The Central Hermite method shows 
the largest errors, while Upwind Hermite shows the smallest errors. The error for the Vir¬ 
tual and Dual methods lie in-between, with the Dual Hermite displaying smaller errors 
at low N. Smaller values of the CFL constant C increase the dispersion and dissipation 
error, though the effect is less noticable as N increases. 

For each method, the error E^h is observed to follow 

_e-ikcdt\ /khV^^^ 

e-ikcdt\ ~ ^ J 

where N is the number of degrees of freedom per node, and the underlying Hermite 
approximation space is of degree 2N+2^ 

'''If we seek instead the error in the discrete and exact wavenumbers |a; —|, we recover convergence rates 
of and for the real and imaginary parts, respectively In comparison, DG with co¬ 
volume filtering achieves rates of 0(/z^^+^) and for the real and imaginary parts, respectively 

while Galerkin methods result in rates of either 2N+1 or 2N+3 (for N even or odd) under a degree N 
approximation space Q. 
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3.2.1 Optimizing dispersive and dissipative errors 

Finally, motivated by Dispersion Relation Preserving (DRP) finite difference schemes 
dispersive and dissipative errors may be improved through optimization of entries 
of the interpolation matrix. As mentioned in Section |2.2[ the Virtual Hermite method 
produces a degree 2N +1 reconstruction using degree N data from three nodes, though 
there is sufficient data to define a higher 3N+2 degree reconstruction. We define the 
interpolation matrix H G ^ (^^^+3) 


H 


H 

H2 


where H is the Virtual Hermite interpolation matrix. DRP schemes enforce a fixed order 
of approximation for a given finite difference stencil, while using additional degrees of 
freedom to optimize the dispersion relation. Similarly, fixing the first 2N+2 rows of H, 
the reconstruction implied by H is enforced to match that of the Virtual Hermite recon¬ 
struction for the first 2N+2 coefficients, while entries of the matrix H 2 (which determine 
higher order coefficients) are used to minimize dispersion and dissipation errors. The 
entries of H depend on the ratio between L^ — x (or R^ —x) and as the grid spacing h. 
Since this ratio is constant as a function of /z, H 2 does not change drastically as a grid is 
refined. However, the optimization does appear to be sensitive to the value of C. 

To demonstrate the effect of optimization, we compare the Virtual Hermite method 
to an optimized scheme for N = 1 and CFL constant C = .9. We produce the optimized 
submatrix H 2 by minimizing the real and imaginary parts of the relative dispersion error 
for the advection equation 

/Re (A/, V /Im (A/, \ ^ 


Re(e 


—ikcdt\ 


+ 


Im(e 


—ikcdt) 


with c = 1 and K = 8 grid cells. The same optimized submatrix H 2 is then used on a finer 
K = 16 grid, and computed solutions for the Virtual and optimized Hermite methods 


in Figure 10 The spectra of the update 


are compared for the initial condition ^-4sm(7rx)2 
matrix S and dispersion/disspation errors are also compared in Figure [TO} The dis¬ 
persion error and spectra are shown to be significantly improved, and numerical results 
indicate that under-resolved features are convected with greater accuracy. 

Unfortunately, the benefits of such an approach appear to be limited to low orders 
of approximation. At higher orders, optimization did not reduce the dispersion and dis¬ 
sipation error significantly compared to the unoptimized scheme. Additionally, the 
stability of such an approach is not guaranteed for Hermite methods (compared to DRP 
schemes, which optimized over symmetric stencils to guarantee stability). For example, 
for N = 1, the spectral radius of the update matrix for the unoptimized scheme was com¬ 
puted to be | 0 (S) = 1 to machine precision. For the N=1 optimized scheme, p{S) 1.0005, 
and strict enforcement of jO(S) < 1 resulted in either non-convergence of the optimization 
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Figure 10: Spectra, dispersion/dissipation errors, and convection of the initial condition 
c“ 4 sin( 7 rx )2 5 periods. Results are shown for both standard and optimized Virtual 
Hermite methods with N = 1, K = 16, and C = .9. Optimization is done on a coarse K = 8 
mesh. 

problem or subpar dispersion and dissipation properties. Further study is required to 
address these issues. 


4 Extension to two dimensions 


Each Hermite method may be extended to higher dimensions naturally through a tensor 
product construction. In this work, we take the grid O to be the tensor product of one¬ 
dimensional grids. Assuming grid spacings hx^hy in the x and y directions, respectively, 
each point {xm,ym) ^ admits the tensor-product expansion 


N N 
j=0k=0 


X-Xn 


f y Vm \ 

\ hy ) ' 


For linear autonomous equations, a temporal Taylor series may be used to evolve the so¬ 
lution in time. Algorithm [^describes this process for the two-dimensional scalar advec- 
tion equation, using derivative matrices Dx,Dy for the x and y coordinates, respectively. 
Due to the tensor-product nature of the Hermite interpolants in higher dimensions, the 
Taylor series must be of order d{2N+l) to be exact in d-dimensions [|15]. Numerical ex¬ 
periments indicate that reducing the degree of the Taylor expansion in time results in a 
tighter timestep restriction; however, this only decreases the restriction by some constant 
factor, which is independent of the order of approximation. In all experiments, the in¬ 
crease in the order of the Taylor expansion did not correspond with a significant decrease 
in error. 

While time evolution is extended in a straightforward way regardless of spatial di¬ 
mension, interpolation operators in higher dimensions are defined through applications 
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Algorithm 2 Time evolution procedure for 2D scalar advection. N may be taken to be 
d{2N+l) for exact time evolution. 

1: procedure TWO-DIMENSIONAL TEMPORAL TAYLOR SERIES EVALUATION 
2: W = u” 

3: for£ = N,]V-l,...,0do 

4: w^-w+^(-cDx-cDy)u" 

5: u"+^=W 


of ID interpolation operators along each coordinate direction. The application of opera¬ 
tors for the Virtual and Central Hermite methods is illustrated in Figure and we refer 
the reader to for more details on the extension of the Dual Hermite method to mul¬ 
tiple dimensions. The Upwind Hermite reconstruction may be adapted to the advection 
equation in higher dimensions by considering the direction of advection along each co¬ 
ordinate, though the direction of the reconstruction will depend on the sign of Cx at each 
point. 



(a) x-reconstmction (Virtual) 
O © O 

o m o 

o @ o 

(d) y-reconstruction (Virtual) 


©-► -© 

o • o 

©--© 

(b) x-reconstruction (Central) 

o © o 

o • o 

o © o 

(e) y-reconstruction (Central) 


O O O 

©-O 

©- ►(§) o 

(c) ^^-reconstruction (Upwind) 

O O O 

O (§) O 

O © O 

(f) y-reconstruction (Upwind) 


Figure 11; Two-dimensional Hermite reconstruction stencils. The Upwind reconstruction 
assumes Cx,Cy > 0. Nodes which contribute data to a reconstruction are circled. 
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4.1 Numerical experiments in two dimensions 


We consider two model problems in two space dimensions; the periodic advection equa¬ 
tion 


du du du 

where c = {cxfCy) is a unit vector, and the isotropic wave equation in first order form 


1 dp 
dt 


—Vu, 


dvL 

dt 


= -Vp, 


where c is a specified wavespeed, p is pressure and u = {u,v) is the velocity. The CFL 
condition for the two-dimensional advection equation is 

\\c\\dt<h, 


while the wave equation depends on the maximum wavespeed in each coordinate di¬ 
rection. For the non-dimensional isotropic wave equation above, this results in the CFL 
condition cdt<h. 

We note that behavior of the Upwind Hermite method is reported only for the advec¬ 
tion equation. While the Upwind Hermite method may be extended to hyperbolic sys¬ 
tems through a characteristic-based approach, numerical experiments with the isotropic 
wave equation indicated that the method resulted in a timestep restriction depending 
on the degree of approximation. We intend to explore additional generalizations of the 
Upwind Hermite method to systems of equations in multiple dimensions in the future. 

As in one dimension, we introduce a CFL constant C such that dt — Ch, and examine 
the behavior of each Hermite method at various values C. Figure[^shows the conver¬ 
gence of the Virtual and Central Hermite methods for the advection equation. We take 
the advection speeds Cx,Cy and exact solution to be 

Cj = cos(7r/3), Cy = sin(7r/3), u{x,y,t)—sirv{n{x — Cxt))sirv{n{y — Cyt)). 

and compute errors for isotropic grids of 8 x 8,16 x 16, and 32 x 32 nodes. 

Unlike Hermite methods in one space dimension, a temporal Taylor series of degree 
d{2N+l) is required for exact time evolution in d dimensions | (T^ . All experiments use 
exact time evolution; however, decreasing the order of the temporal Taylor series to 2N+ 
1 did not result in a significant decrease in error, though the stable timestep restriction 
decreases by a factor of .5. rates of convergence are reported in Table|^for T = 1. As in 
the one-dimensional case, the error is observed to converge at a rate close to 

We repeat the Dual, Virtual, and Central Hermite convergence experiments for the 
periodic wave equation in 2D, using the exact solution 

p{x,y,t) — sin( TTx) sin( zry ) cos {V2nt ). 


Upwind Hermite results are not reported, since a straightforward application of the two- 
dimensional upwind Hermite reconstruction does not yield a stable procedure for the 
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(c)C = .9 


Figure 12: Convergence of errors for each Hermite scheme for the two-dimensional 
advection equation with a smooth sinusoidal solution. 



C = .l 

C = .5 

n 

II 

N 

1 

2 

3 

1 

2 

3 

1 

2 

3 

Dual 

2.90 

5.00 

7.01 

2.97 

5.02 

7.02 

3.05 

5.15 

7.00 

Virtual 

2.83 

5.00 

7.01 

2.95 

5.02 

7.03 

3.04 

5.04 

7.06 

Central 

2.43 

4.96 

7.03 

2.79 

5.00 

7.04 

2.95 

5.08 

7.04 

Upwind 

2.96 

4.98 

6.99 

2.98 

5.02 

7.01 

3.07 

5.14 

7.14 


Table 3: rates of convergence of the Dual, Virtual, and Central Hermite methods for 

the advection equation in two dimensions. 
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wave equation. Convergence experiments are repeated for the set of grids used for ad- 
vection, and Figure plots the errors in p at time T —1 for the Dual, Virtual, and 
Central Hermite methods. The Dual and Virtual Hermite methods produce errors of 
very similar magnitude, while the error for the Central Hermite method is larger by a 
factor of roughly 2^ as observed in ID. Surprisingly, at C = .9 and N = 3, the error for the 
Virtual Hermite method is lower than that of the Dual Hermite method. The rates of 
convergence are reported in Table 




(c)C = .9 


Figure 13: Convergence of errors for the Dual, Virtual, and Central Hermite schemes 
for the isotropic wave equation in two dimensions. 
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C = .l 

C = .5 

n 

II 

N 

1 

2 

3 

1 

2 

3 

1 

2 

3 

Dual 

2.86 

4.93 

6.79 

2.84 

4.73 

6.92 

2.91 

4.77 

7.02 

Virtual 

2.85 

4.93 

6.84 

2.75 

4.96 

7.03 

2.82 

5.07 

6.83 

Central 

2.51 

4.71 

6.05 

2.56 

4.22 

6.82 

3.05 

4.95 

6.84 


Table 4: rates of convergence for the isotropic wave equation in 2D. 
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^m — 1 ^m+l 

(a) Hermite-DG coupling (4-stage RK with 2 
substeps) 


I I EH ® ® ^n+l 



(b) DG-Hermite coupling (patch recovery) 


Figure 14: Coupling between Hermite and DG methods without staggered grids (DG 
nodes are squares, while Hermite nodes are circles). 


4.2 Coupling with Discontinuous Galerkin methods 

Hermite methods may also be coupled to Discontinuous Galerkin methods in order to 
tackle more complicated geometries and boundary conditions. In Q, coupling condi¬ 
tions between DG and the Hermite method are constructed for both the primary and 
auxiliary grids using a least squares reconstruction and high order finite difference sten¬ 
cils. Since the Hermite methods introduced in this work do not require staggered grids, 
the transfer of information is simplified. We will refer to the order of approximation 
for the DG method as m. Hermite methods may transfer information to DG methods 
through the numerical flux. Due to a timestep restriction of 0(/z/m^) for DG compared 
to the 0{h) timestep restriction for Hermite methods, multiple DG substeps must be 
taken for each Hermite timestep. The following numerical experiments use a 4th-order 
Runge-Kutta scheme with 5 stages, and necessitates the evaluation of the numerical flux 
for each stage. To maintain high order convergence, we compute flux contributions by 
evaluating the high order Hermite interpolant in time, as shown in Figure The cou¬ 
pling from DG to Hermite is more fragile, as high order derivative information must be 
determined from the DG solution. In the following experiments, these derivatives are 
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provided via a patch reconstruction at Hermite nodes p^ . Projection onto the Hermite 
basis directly yields high order Hermite coefficients; alternatively, these coefficients may 
then be determined by taking derivatives of the reconstructed polynomial. 

We compute errors and convergence rates using the smooth solution 

=sin(27rx)sin(27ri/)cos(27rt). 

for wavespeed c — ^/l/2. On non-overlapping grids, a polynomial of degree 2N +1 is 
constructed using both Hermite and DG data on neighboring elements. Since the best 
possible convergence rate for DG is take the order of approximation for 

DG to be m = 2N (where N is the degree of the Hermite method) in order to preserve the 
O (2N+1) convergence rate. 




Figure 15: (Left) Coarsest coupled mesh used for coupling Virtual Hermite and DG (Her¬ 
mite nodes are circled). (Right) Convergence of errors in 2D for the isotropic wave 
equation. 

Results are shown in Figure for N = 1,2,3 for the coupled Virtual Hermite-DG 
method, using C = .8 for the Hermite timestep. For N>2, convergence was limited by the 
order of the DG Runge-Kutta scheme, and a smaller timestep must be taken to recover 
optimal convergence rates. Similar observations were made when coupling Dual Her¬ 
mite and DG methods Q. Similar behavior is observed when coupling Central Hermite 
and DG, though the Hermite error increases slightly. 

Unfortunately, the approximation Hermite coefficients using DG becomes less ac¬ 
curate at higher orders, as roughly an order of convergence is lost per derivative with 
patch recovery methods. A salient alternative to patch recovery is Smoothness-Increasing 


^ Optimal convergence rates for upwind DG are typically observed in practice, and are provable on specific 
classes of meshes However, on general meshes, DG methods can expect at most 
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Accuracy-Conserving (SIAC) postprocessing |[7[[^, which produces smooth reconstruc¬ 
tions of the solution which converge with rate Under such a method, opti¬ 

mal convergence rates could be preserved using DG and Hermite methods with degrees 
m=N. The postprocessing of higher order derivatives may also yield additional accuracy 
in Hermite coefficients 


5 Conclusions and future work 

We have presented a generalization of Hermite methods for periodic problems, and have 
investigated two new methods within this framework and compared their performance 
to the original Hermite method in the literature. The original Dual Hermite method re¬ 
sults in a two-node stencil in one space dimension, and requires time integration on both 
primal and staggered (dual) grids. The Virtual Hermite method increases the stencil to 
three nodes, but avoids explicit storage of staggered grid degrees of freedom by fusing 
operations on the auxiliary and primary grid together. The Central Hermite method 
modifies the interpolation procedure in order to avoid a staggered grid, and in doing so, 
maintains a two-node stencil and doubles the timestep restriction. However, to achieve a 
specific error resolution, the Central Hermite method requires almost as many degrees as 
the original Dual Hermite method. Additionally, the Central Hermite method may suffer 
from the propagation of spurious modes. The Upwind Hermite method achieves a res¬ 
olution close to that of the Dual Hermite method, while maintaining the same two-node 
stencil and timestep restriction as the Central Hermite method. However, the stability of 
Upwind Hermite schemes does not appear to generalize in a straightforward manner to 
systems of equations in higher dimensions. Finally, since the Virtual, Central, and Up¬ 
wind Hermite methods do not require dual grids, the coupling between Hermite and DG 
methods is simplified. 

Future work will address variable coefficient problems and explore stable extensions 
of upwind Hermite reconstructions to multi-dimensional wave problems, and to use 
these simplified methods to produce efficient many-core parallel implementations in two 
and three space dimensions. 


6 Acknowledgments 

The authors wish to acknowledge the Matlab codes of the Hermite training library CHIDES 
(http;//www.chides.org), as well as helpful discussions with Daniel Appello. 

Arturo Vargas is supported by an NSF graduate fellowship. Jesse Chan and T. War- 
burton are supported by NSF (award number DMS-1216674). Thomas Hagstrom is sup¬ 
ported by NSF (award number DMS-1418871). 


31 


References 


[1] Mark Ainsworth. Dispersive behaviour of high order finite element schemes for the one¬ 
way wave equation. Journal of Computational Physics, 259:1-10,2014. 

[2] Daniel Appelo, Matthew Inkman, Thomas Hagstrom, and Tim Colonius. Hermite methods 
for aeroacoustics: Recent progress. In 17th AIAA/CEAS Aeroacoustics Conference (32nd AIAA 
Aeroacoustics Conference), Portland, Oregon, 2011. 

[3] Claudio Canuto, M Youssuff Hussaini, Alfio Quarteroni, and Thomas A Zang. Spectral Meth¬ 
ods: Fundamentals in Single Domains. Springer Berlin Heidelberg, 2006. 

[4] Xi Ronald Chen, Daniel Appelo, and Thomas Hagstrom. A hybrid Hermite-discontinuous 
Galerkin method for hyperbolic systems with application to maxwell's equations. Journal of 
Computational Physics, 257:501-520, 2014. 

[5] Prince Chidyagwai, Jean-Christophe Nave, Rodolfo Ruben Rosales, and Benjamin Seibold. 
A comparative study of the efficiency of jet schemes. arXiv preprint arXiv:1104.0542, 2011. 

[6] Bernardo Cockburn, Bo Dong, and Johnny Guzman. Optimal convergence of the original 
DG method for the transport-reaction equation on special meshes. SIAM Journal on Numeri¬ 
cal Analysis, 46(3):1250-1265,2008. 

[7] Bernardo Cockburn, Mitchell Luskin, Chi-Wang Shu, and Endre Siili. Enhanced accuracy by 
post-processing for finite element methods for hyperbolic equations. Mathematics of Compu¬ 
tation, 72(242):577-606, 2003. 

[8] Bernardo Cockburn and Chi-Wang Shu. The local discontinuous Galerkin method for time- 
dependent convection-diffusion systems. SIAM Journal on Numerical Analysis, 35(6):2440- 
2463,1998. 

[9] Richard Courant, Eugene Isaacson, and Mina Rees. On the solution of nonlinear hyperbolic 
differential equations by finite differences. Communications on Pure and Applied Mathematics, 
5(3):243-255,1952. 

[10] Michel O Deville, Paul E Eischer, and Ernest H Mund. High-order methods for incompressible 
fluid flow, volume 9. Cambridge University Press, 2002. 

[11] Evan T Dye. Performance analysis and optimization of Hermite methods on Nvidia GPUs 
using CUD A. Master's thesis. University of New Mexico, 2015. 

[12] Richard S Ealk and Gerard R Richter. Explicit finite element methods for symmetric hyper¬ 
bolic equations. SIAM Journal on Numerical Analysis, 36(3):935-952,1999. 

[13] Karsten Eischer. Convective difference schemes and Hermite interpolation. International 
Journal for Numerical Methods in Engineering, 12(6):931-940,1978. 

[14] Bengt Eornberg. The pseudospectral method: Comparisons with finite differences for the 
elastic wave equation. Geophysics, 52(4):483-501,1987. 

[15] John Goodrich, Thomas Hagstrom, and Jens Lorenz. Hermite methods for hyperbolic initial¬ 
boundary value problems. Mathematics of computation, 75(254):595-630,2006. 

[16] Thomas Hagstrom, Daniel Appelo, Tim Colonius, Matthew Inkman, and Chang Youn Jang. 
Simulation of compressible flows using Hermite methods. The Journal of the Acoustical Society 
of America, 131(4):3429-3429, 2012. 

[17] Jan S Hesthaven and Tim Warburton. Nodal discontinuous Galerkin methods: algorithms, anal¬ 
ysis, and applications, volume 54. Springer, 2007. 

[18] Chang Young Jang, Daniel Appelo, Tim Colonius, Thomas Hagstrom, and Matthew Inkman. 
An analysis of dispersion and dissipation properties of Hermite methods and its application 
to direct numerical simulation of jet noise. In 18th AIAA/CEAS Aeroacoustics Conference (33rd 
AIAA Aeroacoustics Conference), page 2240, 2012. 



32 


[19] Andreas Klockner, Tim Warburton, Jeff Bridge, and Jan S Hesthaven. Nodal discontinuous 
Galerkin methods on graphics processors. Journal of Computational Physics, 228(21):7863- 
7882, 2009. 

[20] P Lesaint and PA Raviart. On a finite element method for solving the neutron transport 
equation. Publications mathematiques et informatiques de Rennes, (S4):l-40,1974. 

[21] Stefano Markidis, Jing Gong, Michael Schliephake, Erwin Laure, Alistair Hart, David 
Henty, Katherine Heisey, and Paul Fischer. OpenACC acceleration of the nekSOOO spec¬ 
tral element code. International Journal of High Performance Computing Applications, page 
1094342015576846, 2015. 

[22] David S Medina, Amik St-Cyr, and Timothy Warburton. High-order finite-differences on 
multi-threaded architectures using OCCA. arXiv preprint arXiv:1410.1387, 2014. 

[23] Jens Markus Melenk and S Sauter. Wavenumber explicit convergence analysis for Galerkin 
discretizations of the helmholtz equation. SIAM Journal on Numerical Analysis, 49(3):1210- 
1243, 2011. 

[24] Jean-Christophe Nave, Rodolfo Ruben Rosales, and Benjamin Seibold. A gradient- 
augmented level set method with an optimally local, coherent advection scheme. Journal 
of Computational Physics, 229(10):3802-3827, 2010. 

[25] Jianxian Qiu and Chi-Wang Shu. On the construction, comparison, and local characteris¬ 
tic decomposition for high-order central WENO schemes. Journal of Computational Physics, 
183(l):187-209,2002. 

[26] Philip J Rasch and David L Williamson. On shape-preserving interpolation and semi- 
Lagrangian transport. SIAM journal on scientific and statistical computing, 11(4):656-687,1990. 

[27] Yu-Xin Ren, Hanxin Zhang, et al. A characteristic-wise hybrid compact-WENO scheme for 
solving hyperbolic conservation laws. Journal of Computational Physics, 192(2):365-386,2003. 

[28] KV Roberts and NO Weiss. Convective difference schemes. Mathematics of Computation, 
pages 272-299,1966. 

[29] Jennifer Ryan, Chi-Wang Shu, and Harold Atkins. Extension of a postprocessing tech¬ 
nique for the discontinuous Galerkin method for hyperbolic equations with application to 
an aeroacoustic problem. SIAM Journal on Scientific Computing, 26(3):821-843, 2005. 

[30] Jennifer K Ryan and Bernardo Cockburn. Local derivative post-processing for the discon¬ 
tinuous Galerkin method. Journal of Computational Physics, 228(23):8642-8664, 2009. 

[31] Benjamin Seibold, Jean-Christophe Nave, and Rodolfo Ruben Rosales. Jet schemes for ad¬ 
vection problems. arXiv preprint arXiv:1101.5374, 2011. 

[32] Christopher KW Tam and Jay C Webb. Dispersion-relation-preserving finite difference 
schemes for computational acoustics. Journal of computational physics, 107(2):262-281,1993. 

[33] Timothy Warburton and Thomas Hagstrom. Taming the CEL number for discontinuous 
Galerkin methods on structured meshes. SIAM Journal on Numerical Analysis, 46(6):3151- 
3180, 2008. 

[34] David L Williamson and Philip J Rasch. Two-dimensional semi-Lagrangian transport with 
shape-preserving interpolation. Monthly Weather Review, 117(1):102-129,1989. 

[35] Olgierd Cecil Zienkiewicz and Jian Zhong Zhu. The superconvergent patch recovery and a 
posteriori error estimates. Part 1: The recovery technique. International Journal for Numerical 
Methods in Engineering, 33(7):1331-1364,1992. 



